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We theoretically study a one-dimensional quasi-periodic Fermi system with topological p- wave 
superfluidity, which can be deduced from a topologically non-trivial tight-binding model on the 
square lattice in a uniform magnetic field and subject to a non-Abelian gauge field. The system 
may be regarded a non-Abelian generalization of the well-known Aubry-Andre-Harper model. We 
investigate its phase diagram as functions of the strength of the quasi-disorder and the amplitude of 
the p- wave order parameter, through a number of numerical investigations, including a multifractal 
analysis. There are four distinct phases separated by three critical lines, i.e., two phases with 
all extended wave-functions (I and IV), a topologically trivial phase (II) with all localized wave- 
functions and a critical phase (III) with all multifractal wave-functions. The phase I is related to 
the phase IV by duality. It also seems to be related to the phase II by duality. Our proposed phase 
diagram may be observable in current cold-atom experiments, in view of simulating non-Abelian 
gauge fields and topological insulators/superfluids with ultracold atoms. 

PACS numbers: 71.23.Ft, 73.43.Nq, 67.85.-d, 


I. INTRODUCTION 

The Aubry-Andre-Harper (AAH) model is a workhorse 
for the study of dynamics of particles in one-dimensional 
(ID) quasi-periodic systems^. Over the past few 
decades, it has been extensively used to theoretically un¬ 
derstand the transport and Anderson localization prop¬ 
erties of these interesting systems, revealing a variety of 
transitions between metallic (extended), critical and in¬ 
sulator (localized) phases^r—. Most recently, the AAH 
model has attracted renewed attentions due to its experi¬ 
mental realization in photonic crystals^ 16 and ultracold 
atom a^ — . It has been found to play a non-trivial role 
to characterize emerging topological states of matter— - — 
and the intriguing phenomenon of quantum many-body 
localizatio n 30 ! 31 . 

The AAH model can be formally derived from the re¬ 
duction of a two-dimensional (2D) quantum Hall system 
to a ID chain22. By taking a Landau gauge for the mag¬ 
netic field with <f> flux quanta per unit cell, the motion of 
an electron in a 2D rectangular lattice with tight-binding 
hopping strengths t x and t y can be described by the Hof- 
stadter Hamiltonian (0^- = 27ri0)22 5 


n 2 D = X 






H.c. 


(1) 


which has a translational symmetry in the ^/-direction. 
As a result, the momentum of the motion along the y- 
axis, k y C [0,27r), is well defined. Defining a reduced 
field operator q via aij = e~ lky ici/\/T^, where L y is 
the number of sites along the y- axis, we deduce from 
PL 2 D the AAH model, 


H = V \t x (cj+i^ + H.c.) + t y Vic\ci 

i 


( 2 ) 


where Vi = 2cos(27 xi(f)-\- k y ). The localization properties 
of the AAH model can be easily understood from the par¬ 
ent Hofstadter Hamiltonian. When ty > t x , the electron 
prefers to hop along the ^/-direction and its wave-function 
thus becomes localized in the x-axis. In contrast, when 
t y < t X: the electron motion will extend over the en¬ 
tire x-axis. Indeed, these two cases are related by the 
so-called Aubry-Andre duality, which can be easily de¬ 
rived by considering two different realizations of gauge 
for the magnetic field-2. On the other hand, the topologi¬ 
cal properties of the AAH model can also be understood 
from the 2D Hofstadter Hamiltonian that underlies the 
topologically non-trivial quantum Hall phenomenon. 

In this work, motivated by the recent proposals 
of simulating non-Abelian gauge fields22r^2 and topo¬ 
logical insulators/superfluid s 22 i 2 Z with ultracold atoms, 
we consider a generalized Hofstadter Hamiltonian, ob¬ 
tained by using a TV-component field operator a ij = 

, • • • ,d-^] T in Eq. (pQ) and by replacing the hop¬ 
ping strengths t x and t y with two SU(TV) matrices T x 
and T y . Taking the same reduction to a ID chain along 
the ^-direction (i.e., using af^ = e~ lky ^ cf^ / with 

p = 1, • • • , TV), we then obtain a generalized non-Abelian 
AAH model, 


U = X [(cJ+itlcCi + H.c.) + Vic\fyCi , 

i 


( 3 ) 


where for simplicity we have assumed that Tj = T y and 

have used c i = [c- 1 ^, • • • , c\ N ^] T . 

To be concrete, we consider the simplest non-Abelian 
case with two-component field operator and the SU(TV = 
2) hopping matrices 


T x — t x <j % iN.&y 


( 4 ) 
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FIG. 1: (Color online). Phase diagram of a quasi-disordered 
Fermi system with a p-wave order parameter A and disorder 
strength t y . The energy is in units of the hopping ampli¬ 
tude t x . Four different phases are separated by the critical 
lines AB, AC, and CD: the phase (I) with spatially extended 
wave-functions, the topologically trivial phase (II) with local¬ 
ized wave-functions, the critical phase (III) with fractal wave- 
functions, and another extended phase (IV). The phases (II) 
and (IV) are dual to the phase (I). Multifractal analysis of 
wave-functions has been made for the 21 points in the critical 
phase (III) or on the critical lines (see Table I). 


and 


Ty — ty & Z , (5) 

where a y and a z are the usual 2 by 2 Pauli matrices. 
In the absence of the quasi-periodic disorder term (i.e, 
t y = 0), the model describes a ID topologically non¬ 
trivial insulator or a spinless p-wave superfimd22ri£I. On 
the other hand, in the limit of A = 0, where T x and 
T y commute with each other, the model reduces to the 
standard A AH model for each component. 

We find that the localization and topological properties 
of the AAH model are profoundly affected by the exis¬ 
tence of a nonzero A. Our main results are summarized 
in the phase diagram Fig. [lj There are four distinct 
phases (I-IV) separated by three critical lines AB, AC 
and CD. For a sufficiently large quasi-disorder strength, 
i.e., above the line AB where t y > t x + A the system 
becomes localized and topologically trivia l 21 ! 22 . Below 
the line AC or its duality line CD, where t y < \t x — A|, 
all the wave-functions of the system are extended. The 
area enclosed by the three separation lines is critical and 
all the wave-functions are multifractal. In this study, 
the phase diagram Fig. [T| is theoretically analyzed by a 
number of numerical approaches, including a multifractal 
analysis for the critical area. 

The rest of the paper is organized as follows. In the 
next section (Sec. II), we explicitly write down the 
Schrodinger equation of our ID quasi-periodic system 
and highlight its particle-hole symmetry and duality. In 


Sec. Ill, we determine the phase diagram by calculat¬ 
ing the inverse participation ratio at a given level of ra¬ 
tional approximations and verifying single-particle wave- 
functions in different phases. In Sec. IV, we introduce 
the multifractal approach and present a scaling analysis 
of the structure of wave-functions in the critical regime. 
In Sec. V, we discuss the behavior of edge modes or Ma- 
jorana fermions in different states. Section VI is devoted 
to conclusions and outlooks. 


II. NON-ABELIAN SU(2) 
AUBRY-ANDRE-HARPER MODEL 


For the non-Abelian AAH model with the hopping ma¬ 
trices in Eqs. and ©, we may write the wave-function 


£[ 


UiA )t + 


| 0 >, 


(6) 


then the Schrodinger equation = e\^>) has the ex¬ 

plicit form, 

tx (^i+i T T tyViUi A (v^+i i) — cui,(7 ) 

A {Ui- |-l Vi— i) t x (Vi-\- 1 T Vi— i) tyViVi — £Vi. (8) 

Recall that Vi = 2cos(27ri</> + k y ). Following the stan¬ 
dard routine for investigating localization properties of 
the AAH model, throughout the paper we set k y = 0 
and take an irrational value <fi = (\/5 — l)/2. This inverse 
of the golden mean may be approached by using the Fi¬ 
bonacci numbers F n : <fr = lim n ^ 00 F n -\/F n , where F n 
is recursively defined by the relation ^n+l — F n + F n - 1, 
with Fq = Fi = 1. Thus, in numerical calculations we 
take the rational approximation 

4>~4>n = D^ L . ( 9 ) 

V n 

To minimize the effect of the periodic boundary condi¬ 
tion, we assume that the length of the system is periodic 
with period F n . 

At the n-th rational approximation, the Schrodinger 
equation then becomes periodic with period F n . Accord¬ 
ing to Bloch’s theorem, the wave-function ip can be char¬ 
acterized by a crystal momentum q x C [— 7r/F n , +7r/F n ) 
and satisfies 'ipi+Fn = e iqxFnr ipi (i = 0,1, • • • ,F n — 1). 
Therefore, if we represent ^ as a vector, 

Ip = [w 0 ,U 0 , Ml, «!,••• ,U Fn - 1 ,V Fn -i] T , (10) 

the Schrodinger equation can be solved by finding the 
eigenvalues and eigenvectors of a 2 F n x 2 F n matrix: 
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where 


and 


Ai = 2 t y cos (27r i(j) n ) 


1 0 
0 -1 


B = 




C = 


t x A 
A t x 


exp ( -iq x F n ). 


( 12 ) 

(13) 


(14) 


The energy spectrum consists of 2 F n bands, each of which 
is a function of the crystal momentum q x . Without loss 
of generality, we consider only one state in each band 
by taking q x = 0. Numerically, we solve the eigen¬ 
value and eigenvector problem of the matrix Eq. ([TT1) 
at a given value of n and then perform a scaling analysis 
with increasing n. The irrational limit is reached when 
we extrapolate our numerical results to the scaling limit 
1/n 0. 


A. Particle-hole symmetry 



FIG. 2: (Color online) Mean inverse participation ratio 
(MIPR) as a function of the quasi-disorder strength at two 
p-wave superfluid order parameter A = 0.5 t x (black squares) 
and A = 1.5 t x (red circles). The dashed line and dot-dashed 
lines show the sharp increase of MIPR at phase boundaries. 
Here, we use n — 15 and F n — 987. 


With our specific choice of SU(2) hopping matrices, 
the non-Abelian AAH model has an interesting particle- 
hole symmetry. That is, the model is invariant under the 
particle-hole transformation: 


In other words, for every particle-like solution (iq,^) of 
the Schrodinger Eqs. © and m with energy E > 0, 
there is always a hole-like solution (v*,u*) with energy 
—E. To highlight this particle-hole symmetry, it is use¬ 
ful to interpret the two components of the system as the 
particle and hole components of a spinless p-wave Fermi 
superfluid, in which the parameter A can be conveniently 
identified as a p-wave order parameter. This interpreta¬ 
tion becomes apparent, if we set = Ci and = cj in 
the non-Abelian AAH model. It leads to a Hamiltonian 
that describes a ID p-wave superfluid in a quasi-periodic 
potential, 

H = E [(*x4+Vi + Acj+iCj + H.c.) + t y Vic\ci , 

(16) 

which shares the same energy spectrum as the AAH 
model. The Anderson localization of the Hamiltonian 
Eq. (fT6l) has been recently analyzed by Lang et 
and DeGottardi et al 21 . It was shown that the system 
becomes localized once t y > t x + A. 

B. Duality 

Our non-Abelian AAH model has an interesting dual¬ 
ity, which is very useful to understand the phase diagram. 


In the p-wave Hamiltonian Eq. m, if we make the fol¬ 
lowing replacement for the field operator, 

Ci —d\ if i is odd (17) 

and 

Ci di if i is even, (18) 

and replace </> by </> + 1/2 in the quasi-periodic potential 
Vi, Eq. © keeps the same form, except that t x and A 
are exchanged. Thus, the system is self-dual if t x = A. 

In the absence of the p-wave order parameter (A = 0), 
the AAH model also has Aubry-Andre duality, which re¬ 
lates a wave-function at A = t y /t x = a to the one at 
A = 1/a by a Fourier transformation^. From the phase 
diagram Fig. [lj a similar duality seems to exist if A ^ 0. 
Unfortunately, we are not able to find a simple transfor¬ 
mation to relate the parameters t x and t y for a given A, 
or the parameters A and t y for a given t x . 


III. INVERSE PARTICIPATION RATIO AND 
PHASE DIAGRAM 

A useful quantity in characterizing phase transitions 
of quasi-periodic systems is the inverse participation ra¬ 
tio (IPR). For a given normalized wave-function -0, it is 
given by JT = JT [uf + vf \, which measures the in¬ 
verse of the number of lattice sites being occupied by 
particles. As we shall see, the non-Abelian AAH model 
considered in this work has a pure energy spectrum so 
that extended, critical and localized wave functions do 
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leads to the proposed phase diagram shown in Fig. [lj 

We find four phases (I-IV) that are separated by three 
critical lines AB , AC and CD. Both phases I and IV have 
extended wave-functions and are related to each other by 
the dual transformation given in Eqs. ( 1171 ) and ([ 151 ) . All 
the wave-functions in the phase II are instead localized. 
It seems that these localized wave-functions may also re¬ 
late to the extended wave-functions in the phase I (or 
IV) by a duality that is analogous to the Aubry-Andre 
duality occurring at A = 0. Yet, such a duality trans¬ 
formation is still to be determined. The three separating 
lines enclose a large area in which all the wave-functions 
are critical. 

In Fig. 03 we examine the representative ground-state 
wave-functions (i.e., of the state at the edge of the en¬ 
ergy spectrum) in different phases with A = 0.5G. With 
increasing strength of the quasi-disorder potential t y /t x , 
it is clear that the wave-function is extended in the phase 
I (see a and b), critical in the phase III (c and d), and 
localized in the phase II (e and f). 


FIG. 3: (Color online) The Bogoliubov quasi-particle wave- 
functions of the lowest edge mode, ui (left panel) and v% (right 
panel), at the p-wave superfluid order parameter A = 0.5 t x . 
From top to bottom, the strength of the quasi-disorder t y 
increases from 0.25 t x (a, b) to t x (c, d), and finally to 1.75 t x 
(e, f), corresponding to the cases of extended, criticial and 
localized wave-functions, respectively. Here, we take n = 15 
and F n — 987. 


not coexist (i.e., there are no mobility edges). There¬ 
fore, at the n-th rational approximation it is convenient 
to define a mean inverse participation ratio (MIPR), 


MIPR = 


1 


2F n -l F n -1 

2 F n ^ ^ 


E 

i E =0 i =0 


(19) 


where ie is the index of energy levels. For extended 
states, MIPR scales like F~ lm , while for localized states, 
MIPR tends to a finite value 0(1). For critical states, 
MIPR behaves like F~ a , where 0 < a < 1 depends 
on the multifractal structure of wave-functions. We use 
MIPR to determine the phase boundaries separating the 
extended, critical and localized phases, which are iden¬ 
tified by the turning points of MIPR as a function of 

ty/t x . 

Figure [2] reports the evolution of MIPR in a logarith¬ 
mic scale at two p-wave order parameters A = 0.5 t x and 
A = 1.5 t x . We have used n — 15 and F n = 987. There 
are two turning points of MIPR located respectively at 
t y = \t x — A| and t x + A, at which MIPR increases very 
rapidly. We have checked that with increasing n, the 
slope of MIPR at these turning points becomes sharper. 
Thus, we anticipate that finally there will be a jump at 
the turning points in the scaling limit n oc, signal¬ 
ing a phase transition. The determination of the turn¬ 
ing points of MIPR at different p-wave order parameters 


IV. MULTIFRACTAL ANALYSIS OF CRITICAL 
WAVE-FUNCTIONS 

To strengthen the proposed phase diagram in Fig. 
we further investigate the scaling behavior of wave- 
functions by using a multifractal analysis^. At the n-th 
level of rational approximations, where the period of the 
lattice is F nj we analyze the probability measure at the 
lattice site i, pi = uj + v\ (i = 0, • • • ,F n — 1), from a 
selected wave-function -0, which is normalized to unity 
JT Pi = 1 • The scaling index on for pi is defined by 

Pi=F~ a '- (20) 

The key observable to characterize the scaling behav¬ 
ior of the wave-function is the singular spectrum f n (a) 
defined by 

Sl n {a) = Fl^ a \ (21) 

where Q n (a)da is the number of lattice sites, which have 
an index cti distributed between a and a+da. The singu¬ 
lar spectrum in the scaling limit f(a) can be calculated 
as /(a) = lim n ^oo / n (a). For extended wave-functions, 
all the lattice sites have a probability measure pi ~ 1 /F n ; 
thus f{a) is only defined at a = 1 with f(a) = 1. For 
localized wave-functions, on the other hand, pi is nonzero 
only on a finite number of lattice sites. These sites have 
an index a = 0 and the remaining sites with exponen¬ 
tially small probability measure have a = oo; thus f(a) 
takes only two values: f(a = 0) = 0 and f(a = oc) = 1. 
For critical wave-functions, the index a has a distribu¬ 
tion and hence the singular spectrum f(a) is a smooth 
function defined on a finite interval [a m in, <a m ax]- There¬ 
fore, it is clear that, to identify extended, critical and 
localized wave-function, we may simply examine the min¬ 
imum value of the index <a, which should take a m [ n = 1 
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FIG. 4: (Color online) Plots of a m in vs 1/n for the wave- 
functions of the lowest edge mode (ie = 0) (a) and of the 
center mode (ie — F n ) (b), at the p -wave superfluid order 
parameter A = 0.5 t x . The black squares, red circles and 
blue triangles show the results with increasing quasi-disorder 
strength t y , as indicated in the left panel (a). 


(extended), 1 > <a m i n > 0 (critical) and <a m i n = 0 (local¬ 
ized), respectively. 

For the numerical calculation of /(a), we follow the 
work by Hiramoto and Kohmotc4 and define an entropy 
function, 


S n (oi) = - lnft n (a), (22) 

n 

which is related to the singular spectrum as (cf. Eq. ([2T]0 
fn (a) = (a), (23) 


where 6 = ln[(\/5 + l)/2]. The entropy function can 
be calculated by an analogous formalism to the usual 
statistical mechanics^. First, one introduces a partition 
function, 

F n ~ 1 

z n(q) = E 

i =0 

and a free energy, 

Qn ( q ) = -\nZ„ ( q ). (24) 

n 

The entropy function is then obtained through the Leg¬ 
endre transformation, 


Sn (a) = Q n (q) + qae 


and 


1 d&n (q) 

e dq 


(25) 


We have calculated f n (a) = S n (ot)/e for finite Fi¬ 
bonacci indices n and have tried to extrapolate them to 


a 

a 


FIG. 5: (Color online) Plots of o m in for the wave-function of 
each mode at A = 0.5 t x and at three different quasi-disorder 
strengths ( t y /t x — 0.25, 1 and 1.75 from top to bottom). We 
have used n — 18 and F n = 4181, so the index of energy mode 
%e runs from 1 to 2 F n — 8362. 



the scaling limit 1/n 0. In Fig. 01 we present exam¬ 

ples of determining a m [ n for wave-functions in the phases 
I, III and II at the p -wave order parameter A = 0.5 t x . 
We consider two typical wave-functions, one at the edge 
of the energy spectrum (denoted by an energy index 
%e — 0) and another at the center of the spectrum 
{ie — F n ). For both wave-functions, a m [ n extrapolates 
to 1 if t y = 0.25 t x , to about 0.58 if t y = t x , and to 0 
if t y = 1.75G- In Fig- El we examine the value of <a m i n 
for all the wave-functions at the n = 18th rational ap¬ 
proximation. For all three disorder strengths, <a m i n varies 
smoothly with increasing spectrum index i#, suggesting 
that there is no mobility edge in the energy spectrum. 
In this way, we confirm that all the wave-functions are 
extended in the phase I, critical in the phase III, and 
localized in the phase II. 

We now focus on the wave-functions in the critical 
phase II and on the three critical lines. Table I summa¬ 
rizes the values of <a m in’s for the edge and the center of 
the energy spectrum for various points in the phase dia¬ 
gram, Fig. [lj At the duality point A of the original A AH 
model (A = 0, t y — t x ), our results of <a m in(edge) ~ 0.17 
and <a m i n (center) 0.36 are consistent with the previous 
calculations by Hiramoto and Kohmoto^. Quite gener¬ 
ally, all the edge states on an individual critical line, ei¬ 
ther AH, AC or CD (excluding the connecting points A 
and (7), seem to have identical values of a m [ n . This is 
also true for the center states, although their values of 
<a m in may be different from those of the edge states. On 
the other hand, all the states in the whole region III, no 
matter at the edge or at the center of the spectrum, have 
identical values of a m [ n 0.58. Therefore, points on an 
individual critical line (excluding the points A and C) or 
in the whole area III may belong to the same universal¬ 
ity class. Further numerical verification of this conjecture 
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TABLE I: Values of a m in for the wave-functions at the edge 
and the center of the energy spectrum. The mark x means 
that there is no enough convergence in data for determining 

C^min • 


A /t x 

ty / tx 

position 

own(edge) 

c^min (center) 

0 

1.0 

A 

0.17 

0.36 

0.25 

1.25 

AB 

X 

0.38 

0.5 

1.5 

AB 

X 

0.39 

0.75 

1.75 

AB 

X 

0.38 

1.0 

2.0 

AB 

X 

0.39 

0.25 

0.75 

AC 

0.58 

0.67 

0.5 

0.5 

AC 

0.56 

0.65 

0.75 

0.25 

AC 

0.59 

0.65 

1.25 

0.25 

CD 

0.65 

0.58 

1.5 

0.5 

CD 

0.64 

0.59 

1.75 

0.75 

CD 

0.66 

0.59 

2.0 

1.0 

CD 

0.64 

0.59 

0.5 

1.0 

III 

0.58 

0.56 

1.0 

0.5 

III 

0.57 

0.57 

1.0 

1.0 

III 

0.57 

0.56 

1.0 

1.5 

III 

0.57 

0.57 

1.5 

1.0 

III 

0.57 

0.58 

1.5 

1.5 

III 

0.58 

0.58 

1.5 

2.0 

III 

0.57 

0.58 

2.0 

1.5 

III 

0.58 

0.57 

2.0 

2.0 

III 

0.57 

0.57 


requires the comparison of the curve f(a) on the whole 
interval [<a m in 5 <a m ax]« Unfortunately, the convergence of 
our numerical estimates for f(a) at arbitrary a is too 
poor to reach a conclusive confirmation. 


V. MAJORANA EDGE MODES 

In the absence of the quasi-periodic potential (t y = 0), 
the Schrodinger equations © and (j8j) hosts zero-energy 
edge mode at the two boundaries, as a result of its non¬ 
trivial topology in the energy band structure. These zero- 
energy edge modes, also referred to as Majorana fermions 
in the case of superfluidity, have been shown to persist at 
nonzero t y: until the localized phase is reache d 21 ! 22 . We 
have calculated the energy of the Majorana edge mode for 
a finite length system, cm, by using the open boundary 
condition, i.e., setting the 2x2 matrix C = 0 in Eq. (HU). 

In Fig. 0 we report the energy gap E gap = 2 cm as 
a function the system length L = F n at A = 0.5G and 
several quasi-disorder strengths. The length dependence 
of the edge-mode energy is dramatically affected by the 
quasi-disorder strength. For small disorder strengths in 
the extended phase I or the critical phase III, the energy 
decreases exponentially as the length of the system L 
increases. On the other hand, on the separating critical 



FIG. 6: (Color online) Length dependence of the energy gap 
at A = 0.5 t x and at different quasi-disorder strengths (as in¬ 
dicated). In the extended and critical phases, the energy gap 
vanishes exponentially with increasing system size, as sug¬ 
gested by the two exponential fitting lines. The inset shows 
the evolution of the energy gap with increasing quasi-disorder 
strength at L — F n — 987, where n — 15. The energy gap 
increases rapidly after entering the localized phase (II). 


line AB or in the localized phase II, the energy decreases 
much slower with increasing length L. At a sufficiently 
large length L = = 987, as shown in the inset of the 

figure, we find that the system opens a nonzero energy 
gap only at the critical line AB , where t y > t x + A, in 
consistent with previous theoretical finding s 21 ! 22 . 


VI. CONCLUSIONS 

In summary, we have proposed a generalization of 
the Aubry-Andre-Harper model to the non-Abelian class 
with SU (N) hopping matrices. The localization and 
topological properties of the model is greatly affected by 
such a generalization. We have performed a systematic 
investigation of the simplest SU(2) case with non-trivial 
p-wave superfluity and have shown that its phase diagram 
becomes much richer. There is a large window for the 
critical phase in the phase diagram, which is separated 
from the extended and localized phases by three criti¬ 
cal lines. Our multifractal analysis of the critical wave- 
functions indicates that points on an individual critical 
line or in the critical phase may belong to the same uni¬ 
versal class. Further issues - such as the spectral statis¬ 
tics on the critical lines and in the critical phase - are of 
interest and will be addressed elsewhere. 

The proposed non-Abelian SU(2) Aubry-Andre- 
Harper model might be realized in cold-atom laborato¬ 
ries in the near future. In particular, in view of recent 
numerous attempts for creating non-Abelian gauge fields 
with ultracold atoms, various non-Abelian Aubry-Andre- 
Harper models could be simulated. We anticipate even 
richer phase diagrams with, for example, non-pure en- 







































7 


ergy spectrum, in which the extended, critical and local¬ 
ized states may coexist and be separated by some mo¬ 
bility edges. Experimentally, it would be interesting to 
observe mobility edges in one-dimensional quasi-periodic 
systems^. 
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